Numerical models of travel-time inhomogeneities in sunspots 
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ABSTRACT 

We investigate the direct contribution of strong, sunspot-like magnetic fields 
to helioseismic wave travel-time shifts via two numerical forward models, a 3D 
ideal MHD solver and MHD ray theory. The simulated data cubes are ana- 
lyzed using the traditional time-distance center-to-annulus measurement tech- 
nique. We also isolate and analyze the direct contribution from purely thermal 
perturbations to the observed travel-time shifts, confirming some existing ideas 
and bring forth new ones: (i) that the observed travel-time shifts in the vicinity 
of sunspots are largely governed by MHD physics, (ii) the travel-time shifts are 
sensitively dependent on frequency and phase-speed filter parameters and the 
background power below the pi ridge, and finally, (iii) despite its seeming limita- 
tions, ray theory succeeds in capturing the essence of the travel-time variations 
as derived from the MHD simulations. 

Subject headings: Sun: helioseismology — Sun: magnetic fields — Sun: oscilla- 
tions — sunspots 
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INTRODUCTION 



Local helioseismic diagnosti c methods such as time -distance helioseismology (IDuva 



1993 ). helioseismic holography ( Lindsey fc Braunlll997 ) and ring-diagram analysis ( Hill 



1 et al. 



19881 ) 



have over the years provided us with unprecedented views of the structures and flows under 
sunspots and active regions. However, a growing body of evidence appears to suggest that 
interpretations of the measured statistical changes in the properties of the wave-field may be 
rendered inaccurate by complexities associated with the observations and wave propagation 
physics. Incorporating the full MHD physics and understanding the contributions of phase 
and frequency filters, and differences in the line formation height, are thought to be central 
to future models of sunspots. 

One of the earliest studie s that highlighted t he interaction of waves with sunspots was 
the Fourier- Hankel analysis of iBraun et al.l (119871 ). who found that sunspots can absorb up 
to half of the incident acoustic-wa ve power and shift the phases of interacting waves quite 
significantly (see also Braun 1995h . These results were echoed over the years by a steady 



steam of theoretica l results (e.g . iBogdan et al.lll996l ; ICally fc Bogdanlll997l ; ICally et al.ll2003 



Crouch et al.ll2005l ; ICally! 120071 ) that have consistently emphasized the need for more sophis- 
ticated modeling and interpretation of wave propagation in strongly magnetized regions. 

Important advanc es i n our observationa l understanding of sunspots were also achieved 



by lDuvall et a 



(19 96T) andlZhao et al.l (120011 ). who inferred the presence of flows underneath 



sunspots, and iKosovichev et al.l (120001 ) who estimated the sub-surface wave-speed topology. 
However, while the inversion procedures applied to derive these results fail to directly account 
for the tensorial nature of magnetic field effects, the action of the field is mimicked via 
changes in the acoustic properties of the medium (the so -called wave speed). Recently 
howe ver, numer ical forward models of h elioseismic wave (e.g. ICameron et al.ll2008l ; iHanasoge 
20081 ) and ray (IMoradi fc Callyl 120081 ) propagation in magnetized atmospheres have been 
developed and a r e beg inning to make inroads into this problem. In particular, the results of 
Moradi fc Callyl (120081 ) and Cameron (2008; private communication) strongly suggest that 
active-region magnetic fields play a substantial role in influencing the wave field, and that 
the complex interaction of magnetic fields with solar oscillations, as opposed changes in the 
wave-speed, are the major causes of observed travel-time inhomogeneities in sunspots. 

We study the impact of strong magnetic fields on wave propagation and the conse- 
quences for time-distance helioseismology using two numerical forward models, a 3D ideal 
MHD solver and MHD ray theory. The simulated data cubes are analyzed using the tra- 
ditional surface-focused ce nter-to- annulus meth od frequently applied in the time-distance 
a nalyses of sunspo t s (e.g ., ICouvidat et al.ll2006l ). Furthermore, we apply the same method 



as 



Moradi fc Callyl (120081 ) to also isolate and analyze the thermal contribution to the observed 
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t ravel-time shifts. 



THE SUNSPOT MODEL 



The background stratification is given by a convectively stable, hydrostatic m = 2.15 
polytrope below the photosphere, smoothly connected to an isothermal atmosphere. The 
magneto-hydrostatic ( MHS) sunspot model t hat we embed in the background is similar in 
construction to that of ICameron et al.l (120081 ) and IHanasogd (120081 ), where the flux tube is 



model ed by an axisymmetric magnetic field geometry based on the ISchliiter fc Temesvary 



( 119581 ) self-similar solution. This approximation requires the following choices for the radial 
(B r ) and vertical (B z ) components of the magnetic field: 



B, 



Mifj{z)e- r2 '^ z) , B r 



-r^(z) 

2 dz 



(i; 



where r and z are the radial and vertical coordinates in cylindrical geometry. The term M 
controls the magnitude of the magnetic field (peak field strength of 3000 G at the photo- 
sphere) and the flux (ttM). The horizontal extent of the tube and its rate of divergence with 
altitude is set by i/)(z) (see Figured] a). 

U pon solving the MHS equations of pressure and Loreritz su pport (described in de- 
tail in ICameron et al.l 120081 : IHanasogd 120081 : iMoradi &: Callyl 120081 ) . we obtain the altered 
thermodynamic stratification of the underlying magnetized plasma (Figure [T] c). 



3. FORWARD MODELS 



3.1. MHD Wave-Field Simulations 



The li neari zed ideal MHD wave equations are integrated according to the recipe of 



Hanasogd (120081 ). We implement periodic horizontal boundaries and place damping sponges 
adjacent to the vertical boundaries to enhance the absorption and transmission of outgoing 
waves. Waves are excited via a pre-computed determini stic source funct i on tha t acts on the 
vertical momentum equation. In a manner similar to iHanasoge et al.l (120081 ) . the forcing 
term is also multiplied by a spatial function that mutes source activity in a circular region of 
10 Mm radius to simulate the suppression of granulation related wave sources in a sunspot. 
Figure [1] shows some of the resulting properties of the simulated vertical (Doppler) velocity 
data extracted from the MHD simulations. 
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3.2. MHD Ray-Path Simulations 



Moradi fc Callyl (120081 ) detail the steps involved in using MHD ray theory to model 
helioseismic ray propagation in magnetized atmospheres. Here, we provide a brief description 
of the magneto-acoustic ray tracing procedure. 

The ray paths are computed in Cartesian geometry in the vertical x-z plane assumed 
to contain both magnetic field lines and ray paths. In this case, we only require the 2D 
dispersion relation with the Alfven wave factored out: 



V 



+ c z )u z K z + a 2 c 2 K 2 k\ + c z N' z k 



t2,2 



or 



2 7>-2\ 2 

a z K K 







(2) 



where K = |k|, c represents the sound speed, a the Alfven speed, N 2 is the squared Brunt- 
Vaisala frequency and u 2 is the square of the isothermal acoustic cut-off frequency. The 
remaining term, ku = B .k (where B is the prescribed magnetic field), represents the 
component of the wavevector k parall el to the magnetic field. The full 3D dispersion relation 
is presented in lMoradi fc Callyl (120081 ). The construction of k is completed by spe cifying the 
governing equations of the r ay paths using the zeroth order eikonal approximation (jWeinberg 
19621 ; iMoradi fc Callyl 120081 ) which are solved using a fourth-order Runge-Kutta numerical 
method. The magneto-acoustic rays stay on the fast-wave dispersion branch at all times. 
It should be noted that neither forward model (i.e., §3.11 13. 2j) accounts for the presence of 
flows. 



4. Travel-Time Measurement Scheme 

The simulated Doppler velocity data-cube of §3.11 (extracted at an observational height 
of 200 km above the photosphere) had dimensions of 200 x 200 Mm 2 x 512 minutes, with 
a spatial resolution of 0.781 Mm and a cadence of 1 minute. For the time-distance calcu- 
lations, we compute cross covariances of oscillation signals at pairs of points on the photo- 
sphe re (source at r-i , recei ver at r 2 ) based on a single-skip center-to-annulus geometry (see 
e.g. ICouvidat et aJbood ). We cross correlate the signal at a central point with signals 



averaged over an annulus of radius A = |ra — rjj around that center. Firstly, we filter 
out the /-mode ridge. Subsequently, standard Gaussian phase-speed in conjunction with 
Gaussian frequency filters centered at 3.5, 4.0 and 5.0 mHz with 0.5 mHz band-widths are 



applied in order to study fre quency dependencies of travel times (e.g. iBraun fc Birchl 12006 



Couvidat &: Rajagurull2007l ). The annular sizes and phase-speed filter parameters used in 
estimating the times shown in Figures [2] and [3] (including the central phase speed (v) and 
full width at half-maximum (FWHM) used) are outlined in Table [H The point-to-annulus 
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cross-covariances are fitted by two Gabor wavelets (e.g.. ICouvidat et al.ll2006l ) to extract the 
required travel times. 

In order to compare theory with simulation, we estimate centre-to-annulus mean time 
shifts, Sr mean , using the MHD ray tracing technique of §3.21 for the same sunspot model (§|2]). 
The single-skip magneto-acoustic rays do not require filtering. Instead, they are propagated 
from the upper turning point of their trajectories, in both the positive and negative x di- 
rections, at a prescribed frequency with horizontal increments of 1 Mm across the sunspot. 
The required range of horizontal skip distances are obtained by altering the shooting angle 
at which the rays are initiated. The skip distances are then binned according to their travel 
path lengths, A, while the travel times are averaged across both the positive and negative 
horizontal directions. For both forward models (i.e., §3.11 13.21) . we only concern ourselves 
with the mean phase time shifts. 



5. RESULTS AND ANALYSES 

5.1. Travel-Time Profiles I. MHD Wave-Field Simulations 

Figure [2] shows maps of 5r meari as well as the frequency filtered azimuthal averages of 
8T mean obtained using time-distance center-to-annulus measurements for the measurement 
geometries indicated in Table [H The 5T mean map for A = 6.2 — 11.2 Mm clearly displays 
positive travel-time shifts, reaching a maximum of around 25 seconds at spot center. A 
similar travel-time shift is observed from the azimuthal average of Sr mean when a frequency 
filter centered at 3.5 mHz is applied to the data. We also observe the magnitude of the 
positive ST mean steadily decrease as we increase the frequency filter to 4.0 mHz, with negative 
5r mean starting to appear in the profile, and by 5.0 mHz the travel times observed inside the 
sunspot are completely negative. For the larger annuli, negative time shifts of increasing 
magnitude are consistently observed as we increase the central frequency of the filter. In 
fact, all Srmean maps for A larger than 8.7 — 14.5 Mm that we measured displayed similar 
Sr mean behavior to the 6.2 — 11.2 and 8.7 — 14.5 Mm bins (albeit with smaller time shifts). 

It is important to take note of both the signs of the travel-time perturbations and 
their apparent frequency dependence. Positive 5r mean have traditionally been interpreted as 
indicative of a region of slower wave propagation in the shallow subsurface layers beneath 
the spot, while negative times are of a wave-speed enhancement. So in essence, the <5r mean 
profiles that we have derived from th e simulation would appe a r to indicate a traditio nal "two- 



layered" wave-speed structure (e.g.. iKosovichev et al.l I200CH : ICouvidat et al.ll2006l ) beneath 



the sunspot. However as can be seen in Figure [H the thermal profile of our model atmosphere 
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is a "one-layer" sunspot model (5c 2 /c 2 < 0) and of the order of ~ —40%. Similarly, changes 
in the sub-surface wave speed, (c 2 + a 2 )/cp — 1 (where c is the unperturbed sound speed), 
lie only in the positives ~ 0-250% (not shown here), with the greatest enhancements seen 
near the surface. The large decrease in the sound speed we observe in our model also raises 
the possibility that current methods of linear inversion may lie beyond their domains of 
applicability. 

We also observed that travel times associated with the smallest measurement geometry 
are most sensitive to the phase-speed filter used, i.e., when the phase-speed parameters were 
adjusted to filter all backgr ound power below the pi ridge, negative ST mean were obtained. 



This behavior was noted by lBraun fc Birch! (120081 ). who determined the causative factor to 



be the background power between the p\ and / ridges. It is unsettling that the sign of the 
time shift may be reversed at will, through small changes in the filter width and center. 



5.2. Travel-Time Profiles II. MHD Ray-Path Simulations 



Figure [3] (frames a-c) show the resultant 5r mean profiles derived from the MHD ray 
tracer for identical measurement geometries as used for the time-distance calculations. The 
similarities between between the ray 5r mean profiles and their time-distance counterparts in 
Figure [2] are striking. Firstly, the ray travel-time perturbation profiles contain predominantly 
negative travel-time shifts for all frequencies, albeit with slightly smaller magnitudes. Sec- 
ondly, a similar frequency dependence of Sr mean is also observed. Generally, high frequency 
rays propagated within the confines of a magnetic field are expected to i) travel faster an d 
ii) propagate longer distances than low frequency rays ( Callyl 120071 ; iMoradi fc Callyi 120081 ) . 
However, one significant difference we can observe in these profiles is the absence of any 
positive travel-time shifts for the A = 3.7 — 8.7 Mm bin. This is significant because the 
exclusively negative ST mean we observe across all geometries not only reflects the one-layered 
wave- and sound-speed profiles below the surface, but also highlights the effects that phase- 
speed filtering can have on time-distance measurements (recall that ray calculations require 
no such filtering). Nonetheless, the overall self-consistency between these results and those 
in §5.11 are very encouraging, despite the 2D nature of the ray calculations. 

Given the fact that ray theory appears to succeed in capturing the essence of the travel- 
time variations as derived from the MHD simulation s, we can isolate t he th ermal component 
of the measured 5r mean using the same approach as IMoradi fe Callyl (120081 ) to ascertain the 
contribution to the travel-time shifts from the underlying thermal structure. To do this, we 
essentially re-calculate the ray paths in the absence of the flux tube while maintaining the 
modified sound-speed profile obtained in §[2J The resulting thermal travel-time perturbations, 
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^ T meani would then be purely a result of thermal (sound-speed) variations along the ray path. 

The resulting <5r^ ean profiles, presented in Figure [3] (frames d-f), surprisingly show that, 
even without the magnetic field, ray theory produces negative travel times - the exception 
being for rays propagated at 5.0 mHz. This indicates that the contribution from the under- 
lying thermal structure is significant enough to modify the upper turning point of the ray 
paths, thus shortening the ray travel times. The appearance of negative travel times for a 
model with a decrease in sound speed would appear to be somewhat counterintuitive, since 
from standard ray theory, one would expect negative time shifts with increases in sound 
speed. The most likely explanation for this phenomenon, is that since both the sound speed 
and density differ from the quiet Sun, consequent changes in the acoustic cut-off frequency 
(c/2H, where H is the density scale height) in the near-surface regions of our model modifies 
the the ray path for waves with frequencies less than 5.0 mHz quite significantly, thereby 
causing negative travel-time shifts. 

However, when comparing with the time perturbations derived from calculations that 
include the magnetic field (i.e. Figure [3] frames a-c), MHD effects appear to be dominant 
contributors to the observed time shifts. This is perhaps most evident for A = 3.7 — 
8.7 Mm (Figure [3] d), where the thermal contribution at spot center appears to make up 
approximately 11% of Sr mean (i.e. Figure [3] a) at 5.0 mHz, 18% at 4.0 mHz and 28% at 
3.5 mHz. For the largest bin, we see a similar contribution at 4.0 and 5.0 mHz, but a much 
greater contribution at 3.5 mHz (45% of <5r mean ). 

6. SUMMARY AND DISCUSSION 

Incorporating the full MHD physics into the various forward models used in local helio- 
seismology is essential for testing inferences made in regions of strong magnetic fields. By 
comparing numerical simulations of MHD wave- field and ray propagation in a model sunspot, 
we find that: i) the observed travel-time shifts in the vicinity of sunspots are strongly de- 
termined by MHD physics, although sub-surface thermal variations also appear to affect ray 
timings by modifying the acoustic cut-off frequency, ii) the time-distance travel-time shifts 
are strongly dependent on frequency, phase speed filter parameters and the background 
power below the p\ ridge, and finally iii) MHD ray theory succeeds in capturing the essence 
of center-to-annulus travel-time variations as derived from the MHD simulations. 

The most unsettling aspect about this analysis is that despite using a background strat- 
ification that differs substantially from Model S (Christensen-Dalsgaard et al. 1996) and a 
flux tube that clearly lacks a penumbra, the time shifts still look remarkably similar (at least 
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qualitatively) to observational time- distance analyses of sunspots. Preliminary tests con- 
ducted with different sunspot models (e.g., different field configurations, peak field strengths 
etc.) have also provided similar results. Given the self-consistency of these results, as derived 
from both forward models, it could imply that we are pushing current techniques of local 
helioseismology to their very limits. It would appear that accurate inferences of the internal 
constitution of sunspots await a clever combination of forward modeling, observations, and 
a further development of techniques of statistical wave-field analysis. 
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Table 1. Annuli Radii and Phase- Speed Parameters 



A (Mm) 


v (km/s) 


FWHM (km/s) 


3.7-8.7 


17.71 


11.94 


6.2-11.2 


21.11 


11.94 


8.7-14.5 


24.36 


11.94 
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Fig. 1. — Some properties of the model atmosphere: a) shows the field configuration of the 
sunspot model, depicted as lines of constant magnetic flux, b) represents the power spectrum 
of the simulated Doppler velocity data-cube with the location of the p- and /-mode ridges, c) 
depicts the near-surface sound speed/thermal profile and d) shows a power map normalized 
to the quiet Sun. 
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Fig. 2. — Center-to-annulus time-distance 5r mean maps (no frequency filtering) and azimuthal 
averages for A= 3.7-8.7 (a,d), 6.2-11.2 (b,e) and 8.7-14.5 Mm (c,f) extracted from the MHD 
wave-field simulations of §3.11 Light solid lines represent Gaussian frequency filtering cen- 
tered 3.5 mHz, dashed lines represent 4.0 mHz and bold solid lines represent 5.0 mHz. 
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Fig. 3. — Center-to-annulus ray Sr mean and <5T^ ean profiles for A= 3.7-8.7 (a,d), 6.2-11.2 
(b,e), and 8.7-14.5 Mm (c,f) computed using the MHD ray calculation recipe of §3.21 Light 
solid lines represent 3.5 mHz, dashed lines represent 4.0 mHz and bold solid lines represent 
5.0 mHz. 



